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ABSTRACT 

A 1 -dimensional material model was developed for simulating the transverse 
(thickness-direction) loading and unloading response of aluminum honeycomb 
structure. The model was implemented as a user-defined material subroutine 
(UMAT) in the commercial finite element analysis code, ABAQUS®/Standard. 
The UMAT has been applied to analyses for simulating quasi-static indentation 
tests on aluminum honeycomb-based sandwich plates. Comparison of analysis 
results with data from these experiments shows overall good agreement. 
Specifically, analyses of quasi-static indentation tests yielded accurate global 
specimen responses. Predicted residual indentation was also in reasonable 
agreement with measured values. Overall, this simple model does not involve a 
significant computational burden, which makes it more tractable to simulate other 
damage mechanisms in the same analysis. 


INTRODUCTION 

Honeycomb structure consists of an array of repeating unit cells with a 
geometry such as the hexagonal shape illustrated in Figure 1 [1]. Cell wall 
thickness and cell size are typically varied, yielding honeycombs with a range of 
densities and out-of-plane properties aimed at serving a particular need. 
Honeycombs are manufactured out of a variety of materials, with aluminum and an 
aramid fiber paper, known as Nomex® 3 [2], being the most commonly utilized 
materials in aerospace applications. The relatively low density of honeycomb, 
combined with its high specific out-of-plane compression and shear properties, 
makes this type of structure well suited for weight-critical applications. For 


3 Nomex 1 ' is a registered trademark of E.I DuPont de Nemours, Wilmington, DE, USA. 



example, sandwich structure consisting of honeycomb panels reinforced with stiff, 
thin facesheets, are often used in structural components of aerospace vehicles. 

A drawback, however, of using honeycombs as core materials in sandwich 
structure is the propensity of the honeycomb to crush during a transverse loading 
event (such as low-velocity impact) [3-6]. This damage can result in a significant 
reduction of the in-plane compressive strength of the sandwich panel [3], The 
problem is further exacerbated by the difficulty in detecting such damage, 
particularly after low-velocity impact events that may be barely visible. 
Consequently, some effort has been applied towards experimentally characterizing 
damage that arises from a low-velocity impact event and predicting the resulting 
residual in-plane compressive strength of the damaged panel [7-20]. In many of 
these examples, the simulations used to predict residual compressive strength 
involve an assumed initial damage state in the core material [10, 12, 14, 18] or 
completely disregard the effect of the core [11]. 

Researchers have also focused on predicting core damage in sandwich 
structures arising from localized, transverse loading [20-25]. Many of these efforts 
utilize beam theory solutions whereby a Winkler foundation is employed to 
represent the elastic response of the core material and plastic foundations to 
represent core crushing [21-22], More recently [23-24], a modified beam theory 
analysis was conducted in which a beam is sectioned into segments allowing for the 
progressive development of indentation in the sandwich beam as transverse loading 
is applied. Finite element simulations have also been conducted for simulating low- 
velocity impact on sandwich structures [25]. In this example, dynamic analyses 
were performed, utilizing an elastic-perfectly plastic material parameter to represent 
thickness-direction loading of an aluminum honeycomb core. Other finite element 
simulations of honeycomb core-based sandwich structure [20] have been conducted 
where the cellular geometry of the core was explicitly modeled, and a plasticity 
model was utilized to capture yielding involved during the crush and post-crush 
response of the aluminum honeycomb. These models were applied to simulating 
the quasi-static indentation (QSI) on sandwich plates. 

The vast majority of these simulations utilize an empirically based constitutive 
model for representing the response of honeycomb to transverse loading. These 
models are themselves based on data obtained from a flatwise compression/tension 
test (such as ASTM C365 [26]) of the honeycomb material in question. An 
adequate representation of this response is required in order to predict the extent of 
damage that arises in the honeycomb structure after a localized transverse loading 
event. 

The objective of the work detailed in this paper was to develop a simple one- 
dimensional model for simulating the response of a honeycomb structure to 
transverse loading. As in previous cases, the model is based on the aforementioned 
constitutive relationship. However, the current effort differs from previous work 
because the model is implemented as a user material model into a commercial finite 
element analysis code (ABAQUS® 4 /Standard). Also, the parameters of the 
constitutive model can be easily modified by changing the user input of the material 
model. Finally, this simple approach does not involve a significant computational 
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burden, which makes it more tractable to simulate other damage mechanisms in the 
same analysis. 

The remainder of this paper describes the development of the material model, 
its implementation into ABAQUS®/Standard, and the exercises undertaken to 
evaluate the material model. 


CORE DAMAGE MATERIAL MODEL 

The purpose of the model described in this section is to represent the response 
of aluminum honeycomb structure to transverse loading. The model is based on a 
transverse stress-strain relationship measured using a flatwise compression/tension 
test. This stress-strain relationship is implemented into ABAQUS®/Standard as a 
user-defined, 1 -dimensional material model. Details of the model, its 
implementation into ABAQUS®/Standard and evaluation of the implemented 
model follow in the remainder of this section. 


Constitutive Law 

The stress-strain relationship on which the current honeycomb material model is 
based is typically measured using the flatwise-compression specimen depicted in 
Figure 2. The specimen is based on that detailed in the standard test methods for 
conducting flatwise compression and tension tests (ASTM C365 [26] and ASTM 
C297 [27], respectively), and consists of a 50mm-square block of honeycomb 
bonded between two loading blocks. A typical stress-strain response from a 
flatwise compression test is illustrated in Figure 3. Compressive loading is depicted 
by the solid grey lines in Figure 3, while the solid orange lines indicate unloading. 
Specimens are loaded in compression until the honeycomb cell walls begin to 
collapse, which generally corresponds to the sudden load drop shown in Figure 3. 
The honeycomb specimen continues crushing upon further compressive loading, 
which is characterized either by a stress plateau or a slight increase in stress as 
indicated in Figure 3. The specimen is unloaded after the required amount of core 
crushing has been obtained. This initially results in an elastic response from the 
specimen as the collapsed cell walls begin to defonn. At some point during 
unloading, the specimen undergoes an overall reduction in stiffness, which is 
attributed to plastic deformation in aluminum honeycombs [23]. Complete 
unloading of the specimen either involves fracturing of the cell walls [23], or 
further plastic deformation, resulting in a residual tensile stress [20] as depicted in 
Figure 3. This constitutive relationship is specific to aluminum honeycombs and so 
may well adopt a different form in honeycombs made from other materials. 

The current model idealizes the stress-strain relationship in Figure 3 as a series 
of linear functions depicted by the dashed lines in the figure. The parameters 
associated with this idealization are also depicted in the figure and act as the user 
input for the implementation of this model into ABAQUS®/Standard. 



Implementation of Model into ABAQUS®/Standard 


The constitutive model described in the previous section and illustrated in 
Figure 3 is implemented into ABAQUS®/Standard as a user material model 
(UMAT) [28]. The implemented model can be used with 2-node truss elements 
(ABAQUS®/Standard element type T3D2) in a finite element analysis in which the 
truss elements are used to represent a honeycomb structure. The UMAT is in the 
form of a user-defined subroutine, written in the computer language, 
FORTRAN 77. ABAQUS®/Standard provides the user with access to a range of 
its own subroutines that may be called within the UMAT, enabling access to 
information such as analysis variables (state, field, etc.) and analysis time (step and 
increment number, time step, etc.). The parameters that define the constitutive 
model (Figure 3) are provided by the user in the analysis job input file (.inp) that 
contains all other user-defined information (node numbering, mesh connectivity, 
material, properties, load and boundary conditions, etc.) necessary for running an 
analysis. The UMAT is called at the time of executing an analysis at which point 
the UMAT is compiled and configured for operation within the analysis job. 
Operation of the UMAT within a finite element analysis proceeds as follows: 

1 . Analysis begins and user-input (Figure 3) is read by the UMAT. 

2. The UMAT computes the tangent moduli, Ei, E 2 , and E 3 , (Figure 3). 

3. At the beginning of each analysis increment, the current strain and strain 
increment is provided to the UMAT by ABAQUS®. The UMAT checks 
these values against the strain values that define key stages of deformation 
in the constitutive model (ei, e 2 , and e 3 in Figure 3). The strain increment is 
used to determine the direction of loading in an element. 

4. The UMAT checks the current compressive strain value against a previously 
stored maximum value. If the current value exceeds the stored value then 
this strain is replaced by the current strain. For instance, if the current strain 
is -lOlOpts and the stored value is -1000p,e then the stored value is replaced 
by the value of the current strain. This information is used to store the 
current damage state in an element. 

5. A tangent modulus is selected based on the current strain and strain 
increment. This modulus is assigned to the stiffness matrix of the element. 
The stress in the element is then computed using the selected tangent 
modulus and the current strain increment. 

6. The tangent stiffness and stress values are reported to ABAQUS®. 

7. Steps 5 and 6 are repeated until the analysis converges corresponding to the 
end of the increment. 

8. Steps 3-7 are repeated until completion of an analysis step(s). 

The above operation is perfonned for each truss finite element. Additionally, a 
state of damage is defined in each element that corresponds to one of the three 
damage states depicted in Figure 3 (numbers in parenthesis in Figure 3). Damage 
state 1 corresponds to the region of the constitutive model before cell wall collapse 
has occurred (i.e. the initial elastic response). Damage state 2 corresponds to the 
region of the constitutive model during initial cell wall collapse. Damage state 3 
corresponds to the region of the constitutive model during continued crushing of the 



honeycomb. This information is used by the UMAT to determine the correct 
tangent stiffness during unloading (as illustrated by the dashed arrows in Figure 3) 
and also to provide a record of the damage state of each element at any moment 
during an analysis. 


Evaluation of Material Model 

SIMULATION OF QSI TESTS 

The UMAT was evaluated by applying it to finite element analyses of a 
sandwich plate subjected to quasi-static indentation (QSI) loading. This test case 
was chosen because the specimen response does not involve flexure, which would 
not be captured by the current one-dimensional model. The sandwich plate 
consisted of an aluminum honeycomb core reinforced with quasi-isotropic 
graphite/epoxy composite facesheets. Complete details of the sandwich plate are 
given in Figure 4. QSI tests were previously conducted on this sandwich plate 
using a 25mm-diameter indenter and also a 76mm-diameter indenter [19]. 
Specimens were supported by a rigid back plate during each QSI test. A schematic 
of the test configuration is shown in Figure 4. The indentation/force response of 
each specimen was recorded in addition to the residual indentation in the sandwich 
plate at the end of each test [19]. Finite element analyses of these two test cases 
were conducted as part of an exercise to evaluate the developed UMAT. The 
76mm indenter test case was deemed appropriate for evaluating the UMAT because 
the QSI response of the specimens was found to be largely attributed to core 
damage rather than facesheet damage [19], which is not accounted for in the 
following analyses. Even though significant facesheet damage (in the fonn of 
delamination) was observed in the 25mm indenter case, an analysis of this case was 
conducted nonetheless. In both cases, static, geometrically nonlinear analyses were 
conducted in order to account for indenter contact (as described below). 

Finite element models used in the current evaluation of the UMAT consisted of 
4-node shell elements (ABAQUS® type S4) to represent the upper facesheet. Each 
node of the shell elements was connected to 2-node truss elements, (ABAQUS® 
type T3D2) used to represent the honeycomb core. The length of the truss elements 
was equal to the honeycomb core height (in this case 25mm). A typical finite 
element mesh is shown in Figure 5 (only small number of truss elements are 
illustrated in the figure for clarity). Two-axis symmetry was assumed in the 
analysis in order to reduce the overall model size. Therefore, the center of the 
sandwich specimen corresponds to the comer of the mesh labeled ‘A’ in Figure 5. 
A relatively fine mesh was used to represent the facesheet with an element length of 
0.5mm. All shell elements representing the facesheet had 1:1 aspect ratios. A 
similar meshing scheme was shown to yield a converged solution in a similar 
analysis conducted previously [18], and therefore a mesh convergence study was 
deemed to be unnecessary in the current case. The facesheet that made contact with 
the rigid back plate in actual tests (as illustrated in Figure 4) was initially modeled 
with shell elements. However, this was found to be unnecessary and instead the 
boundary condition detailed in Figure 5 (w=0) was used to represent the supported 
specimen side. The stacking sequence of the upper facesheet was represented using 
a composite stack feature provided by ABAQUS®/Standard, where the ply 



orientation corresponded to the coordinate system shown in Figure 5. The nine 
engineering constants used to define the elastic response of the facesheet plies 
(facesheet damage was not considered in the current analyses) are presented in 
TABLE I, which also includes user input for the UMAT. The UMAT input data 
were obtained from previous flatwise compression/tension tests conducted on the 
sandwich structure in question [20]. The spherical indenters were represented using 
rigid elements (ABAQUS®/Standard type R3D3 and R3D4) and a frictionless 
contact algorithm provided by ABAQUS® was employed to facilitate contact 
between the rigid indenter elements and the shell elements representing the upper 
facesheet. The indenter mesh was positioned such that the node corresponding to 
the indenter tip was coincident with the plane of the facesheet, and the in-plane 
position of the indenter tip corresponded to the center of the sandwich specimen. 
Loading was applied to the rigid indenter mesh by prescribing a displacement along 
the global z-axis to a reference node positioned in the center of the indenter as 
shown in Figure 5. The nodes of the rigid indenter mesh were kinematically 
coupled with this dummy node such that they mimic any translation and rotation 
prescribed at this node. All prescribed boundary conditions are illustrated in Figure 
5. Two analyses were performed for simulating QSI tests conducted using a 25mm- 
diameter and 76mm-diameter indenter, respectively. Each analysis was conducted 
in two steps. The first step involved translating the indenter in the z-axis towards 
the specimen up to a maximum prescribed indentation (denoted as z max in Figure 5). 
In the second step, the prescribed indenter translation was reversed until the 
indenter returned to its position. The maximum allowed time over which a solution 
is sought during an increment (referred to as ‘time increment’ in ABAQUS® 
tenninology), Ati nc _ max , was limited to 0.1% of the total time step in each analysis 
(this limit on Ati nc - m ax was imposed as a result of the study described at the end of 
this section). The computed indentation/force response and residual indentation 
profiles from both analyses were compared with corresponding data measured 
previously [19]. 

MAXIMUM ANALYSIS TIME INCREMENT 

A static, geometrically nonlinear analysis in ABAQUS® is comprised of a 
series of steps such as those detailed above for translating the rigid indenter. In 
turn, each step is partitioned into a number of increments for which ABAQUS® 
attempts to converge to a solution. Typically, the amount of time step over which 
ABAQUS® attempts to converge to a solution, otherwise known as the time 
increment, is set to a relatively small value at the beginning of an analysis step. As 
the analysis proceeds, and if convergence is easily found, ABAQUS® 
automatically increases the time increment up to a maximum time increment that is 
either set by the user or ABAQUS® (see [29] for further details on analysis 
solution strategies utilized by ABAQUS®). It was found that the solution of an 
analysis utilizing the currently developed UMAT was very sensitive to this 
maximum time increment, At; nc . max . This sensitivity was investigated using the 
analysis of the QSI test involving the 76mm-diameter indenter. During this study, 
an analysis of this test configuration was repeated five times, each using a different 
value of Ati nc _ max , ranging from 2% of the total time step to 0.1% of the total time 
step (in each case, the initial time increment was kept constant at 0.1% of the total 



time step). The indentation/force responses and residual indentation profiles were 
computed and checked for convergence relative to Ati nc _ max . 


RESULTS / DISCUSSION 
Simulation of QSI Tests 

The compressive force/indentation responses computed from simulations of the 
QSI tests with a 25mm-diameter indenter and 76mm-diameter indenter are 
presented in Figure 6. The indentation/force response from the actual tests is 
superimposed onto each plot. The overall force/indentation response is captured 
very well by comparison to the measured responses. The beginning of each 
response is linear, corresponding to the initial linear response of the honeycomb. 
After some point, the honeycomb cell walls begin collapsing, which is 
characterized by an initial reduction in overall stiffness of the specimen (stiffness 
here is defined as the gradient to the tangent of the indentation/force curve). The 
amount of crushed honeycomb increases as indentation is increased. However, this 
is accompanied by an increase in contact area of the indenter with the neighboring 
facesheet. Consequently, continued indentation results in a gradual increase in 
specimen stiffness. After the maximum indentation is reached the indenter is 
returned to its original position. The resulting unloading is nonlinear largely 
because the amount of contact area between the indenter and specimen reduces 
upon further unloading. 

The computed residual indentation profiles (along Segment AB in Figure 5) 
from analyses of the 25mm-diameter and 76mm-diameter indenter cases are 
presented in Figure 7. The computed dent profiles are mirrored about the yz plane 
(Figure 5) in order to provide a meaningful comparison with the measured dent 
profiles. The measured residual indentation profiles [19] are superimposed onto 
both plots for comparison. The results show that the analysis of the 25mm indenter 
case overestimated the residual indentation area although the correct overall shape 
of the dent, including the maximum indentation, was reproduced. The computed 
residual indentation from analysis of the 76mm indenter case was very close to the 
measured dent profile. 

The peak applied force, residual indentation and total energy dissipation (area 
encompassed by the force/indentation responses) computed from both analyses are 
presented in TABLE II. The corresponding measured values are also included in 
the table for comparison. In the 25mm-diameter indenter case, the computed peak 
force is almost identical to the measured value (less than 1% difference). The 
computed residual indentation for this case, however, underestimates the measured 
value by almost 19%. The computed total energy dissipation compares more 
favorably and underestimated the measured value by 10%. This significant 
difference in dent profiles is likely attributed to the fact that significant 
delamination was observed in the 25mm indenter case, which was not accounted for 
in the current analysis. 

The peak force computed from the analysis of the 76mm-diameter indenter case 
was again in excellent agreement with the measured value (less than 1% 
difference). The computed residual indentation was in better agreement with the 



measured value (compared to the 25mm-diameter indenter case) and in this case 
underestimated the measured residual indentation by only 6%. The computed total 
energy dissipation underestimated the measured value by just under 10%. This 
relatively improved agreement in the 76mm indenter case is consistent with the fact 
that the specimen response was largely attributed to core damage [19]. 

In general, the overall computed specimen response compares very favorably 
with the experimental data in the 76mm indenter case, and suggests that the 
implemented one-dimensional material model is successful in capturing the 
response of honeycomb structure to transverse loading. 


Maximum Analysis Time Increment 

The maximum allowed time step during each increment, Ati nc - m ax, in the above 
two analyses was limited to 0.1% of total time step. As mentioned previously, this 
limit was imposed as a result of a study conducted to determine the effect of 
maximum allowed time step on analysis results. The results of this study (based on 
analyses of the 76mm-diameter indenter case) are presented in Figure 8, where the 
computed peak force (Figure 8a) and maximum residual indentation (Figure 8b) are 
plotted as functions of the number of increments per analysis step. The number of 
increments per step is the reciprocal of Ati nc _ max , and the results are plotted in this 
manner to highlight the convergence of the solutions. Included in the plots are the 
measured peak force and residual indentation for this test case (denoted by dashed 
horizontal lines in each plot). Both plots show that the computed values tend to 
converge towards their measured counterparts as the number of increments per 
analysis step increases. Hence, from these results it was assumed that a converged 
solution was obtained when the number of increments per step was greater than or 
equal to 1000, i.e. when Ati nc . max 0.001. 

SUMMARY AND CONCLUDING REMARKS 

A material model for representing transverse loading of honeycomb structure was 
implemented as a material model into the commercial finite element analysis code, 
ABAQUS®/Standard. The model is based on an empirically determined stress- 
strain relationship measured using a flatwise/compression specimen. The exercises 
conducted for evaluating the UMAT resulted in the following general observations: 

• The results from analyses of the quasi-static indentation (QSI) tests showed 
overall good agreement with the corresponding test data. This favorable 
comparison indicates that the one-dimensional material model can 
adequately capture the honeycomb core response, at least for QSI 
simulations. 

• Results from analyses that employ the implemented material model would 
be sensitive to the time step taken during each analysis increment. 
However, this study indicated that a converged solution could be found 
when the maximum time increment in an analysis step was limited to 0.1% 
of the total time step. 



In general, this simple approach does not involve a significant 
computational burden, which should make it more tractable to simulate 
other damage mechanisms in the same analysis 
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Figure 1. Hexagonal cell geometry of a honeycomb structure. 



Figure 2. Schematic of a 50mm-square flatwise compression specimen. 
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Figure 3. Idealized stress-strain relationship of honeycomb structure. 



Figure 4. Details of quasi-static indentation test setup [20], 



BOUNDARY CONDITIONS: 

Nodes along Segment AB: 

u=0; <))y=0* 

Nodes along Segment AC: 

v=0; (J) x =0* 

Truss nodes (side opposite facesheet): 

w=0 


Reference node: 

Fixed except for w=z max 

z max =1 -44mm (25mm indenter case) 

z max =2.57mm (76mm indenter case) 

* (fr denotes rotation about the i th axis 


Figure 5. Finite element mesh and boundary conditions of QSI analyses (76mm indentor shown). 
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Figure 6. Computed versus measure force/indentation response of QSI specimens. 
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Figure 7. Computed versus measured residual indentation depths in QSI specimens. 
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Figure 8. (a) Computed peak force versus increments per step, (b) Computed maximum resicual dent 

versus increments per step. 


TABLE I. MATERIAL PROPERTIES USED IN ANALYSES. 


Facesheet ply properties [20] 

El 1 = 143.0 GPa 

E22 = 12.9 GPa 

E 33 = 11.7 GPa 

vi 2 = 0.32 

vi 3 = 0.32 

V23 = 0.44 

Gi 2 = 4.1 GPa 

G 13 =4.1 GPa 

G23 = 3.98 GPa 


UMAT input (Figure 3) [20] 


(T|,=2MPa 

a c =0.84MPa 

a cc =1.12MPa 

a tt =0.01MPa 

a t =lMPa 

k=1.26 

d=l 180p,8 

e 2 =l 0660^,8 

63=98800^8 


TABLE II. COMPUTED AND MEASURED QSI SPECIMEN RESPONSE. 


25mm indentor QSI test 

Computed 

Measured [20] 

% difference 

Peak force (kN) 

1310 

1300 

< 1 % 

Max residual dent (mm) 

0.48 

0.59 

-19% 

Total energy dissipation (Nmm) 

632 

701 

-9.8% 

76mm indentor QSI test 

Computed 

Measured [20] 

% difference 

Peak force (kN) 

2790 

2800 

< 1 % 

Max residual dent (mm) 

1.09 

1.16 

- 6 % 

Total energy dissipation (Nmm) 

2340 

2599 

-9.9% 




